home *** CD-ROM | disk | FTP | other *** search
/ Aminet 30 / Aminet 30 (1999)(Schatztruhe)[!][Apr 1999].iso / Aminet / gfx / misc / gnuplot-3.7src.lha / gnuplot-3.7src / gnuplot-3.7.lha / gnuplot-3.7 / specfun.c < prev    next >
C/C++ Source or Header  |  1998-12-10  |  21KB  |  871 lines
  1. #ifndef lint
  2. static char *RCSid = "$Id: specfun.c,v 1.21 1998/04/14 00:16:20 drd Exp $";
  3. #endif
  4.  
  5.  
  6. /* GNUPLOT - specfun.c */
  7.  
  8. /*[
  9.  * Copyright 1986 - 1993, 1998   Thomas Williams, Colin Kelley
  10.  *
  11.  * Permission to use, copy, and distribute this software and its
  12.  * documentation for any purpose with or without fee is hereby granted,
  13.  * provided that the above copyright notice appear in all copies and
  14.  * that both that copyright notice and this permission notice appear
  15.  * in supporting documentation.
  16.  *
  17.  * Permission to modify the software is granted, but not the right to
  18.  * distribute the complete modified source code.  Modifications are to
  19.  * be distributed as patches to the released version.  Permission to
  20.  * distribute binaries produced by compiling modified sources is granted,
  21.  * provided you
  22.  *   1. distribute the corresponding source modifications from the
  23.  *    released version in the form of a patch file along with the binaries,
  24.  *   2. add special version identification to distinguish your version
  25.  *    in addition to the base release version number,
  26.  *   3. provide your name and address as the primary contact for the
  27.  *    support of your modified version, and
  28.  *   4. retain our contact information in regard to use of the base
  29.  *    software.
  30.  * Permission to distribute the released version of the source code along
  31.  * with corresponding source modifications in the form of a patch file is
  32.  * granted with same provisions 2 through 4 for binary distributions.
  33.  *
  34.  * This software is provided "as is" without express or implied warranty
  35.  * to the extent permitted by applicable law.
  36. ]*/
  37.  
  38. /*
  39.  * AUTHORS
  40.  *
  41.  *   Original Software:
  42.  *   Jos van der Woude, jvdwoude@hut.nl
  43.  *
  44.  */
  45.  
  46. #include "plot.h"
  47. #include "fnproto.h"
  48.  
  49.  
  50. extern struct value stack[STACK_DEPTH];
  51. extern int s_p;
  52. extern double zero;
  53.  
  54. #define ITMAX   200
  55.  
  56. #ifdef FLT_EPSILON
  57. # define MACHEPS FLT_EPSILON    /* 1.0E-08 */
  58. #else
  59. # define MACHEPS 1.0E-08
  60. #endif
  61.  
  62. /* AS239 value, e^-88 = 2^-127 */
  63. #define MINEXP  -88.0
  64.  
  65. #ifdef FLT_MAX
  66. # define OFLOW   FLT_MAX        /* 1.0E+37 */
  67. #else
  68. # define OFLOW   1.0E+37
  69. #endif
  70.  
  71. /* AS239 value for igamma(a,x>=XBIG) = 1.0 */
  72. #define XBIG    1.0E+08
  73.  
  74. /*
  75.  * Mathematical constants
  76.  */
  77. #define LNPI 1.14472988584940016
  78. #define LNSQRT2PI 0.9189385332046727
  79. #ifdef PI
  80. # undef PI
  81. #endif
  82. #define PI 3.14159265358979323846
  83. #define PNT68 0.6796875
  84. #define SQRT_TWO 1.41421356237309504880168872420969809    /* JG */
  85.  
  86. /* Prefer lgamma */
  87. #ifndef GAMMA
  88. # ifdef HAVE_LGAMMA
  89. #  define GAMMA(x) lgamma (x)
  90. # elif defined(HAVE_GAMMA)
  91. #  define GAMMA(x) gamma (x)
  92. # else
  93. #  undef GAMMA
  94. # endif
  95. #endif
  96.  
  97. #ifndef GAMMA
  98. int signgam = 0;
  99. #else
  100. extern int signgam;        /* this is not always declared in math.h */
  101. #endif
  102.  
  103. /* Global variables, not visible outside this file */
  104. static long Xm1 = 2147483563L;
  105. static long Xm2 = 2147483399L;
  106. static long Xa1 = 40014L;
  107. static long Xa2 = 40692L;
  108.  
  109. /* Local function declarations, not visible outside this file */
  110. static double confrac __PROTO((double a, double b, double x));
  111. static double ibeta __PROTO((double a, double b, double x));
  112. static double igamma __PROTO((double a, double x));
  113. static double ranf __PROTO((double init));
  114. static double inverse_normal_func __PROTO((double p));
  115. static double inverse_error_func __PROTO((double p));
  116.  
  117. #ifndef GAMMA
  118. /* Provide GAMMA function for those who do not already have one */
  119. static double lngamma __PROTO((double z));
  120. static double lgamneg __PROTO((double z));
  121. static double lgampos __PROTO((double z));
  122.  
  123. /**
  124.  * from statlib, Thu Jan 23 15:02:27 EST 1992 ***
  125.  *
  126.  * This file contains two algorithms for the logarithm of the gamma function.
  127.  * Algorithm AS 245 is the faster (but longer) and gives an accuracy of about
  128.  * 10-12 significant decimal digits except for small regions around X = 1 and
  129.  * X = 2, where the function goes to zero.
  130.  * The second algorithm is not part of the AS algorithms.   It is slower but
  131.  * gives 14 or more significant decimal digits accuracy, except around X = 1
  132.  * and X = 2.   The Lanczos series from which this algorithm is derived is
  133.  * interesting in that it is a convergent series approximation for the gamma
  134.  * function, whereas the familiar series due to De Moivre (and usually wrongly
  135.  * called Stirling's approximation) is only an asymptotic approximation, as
  136.  * is the true and preferable approximation due to Stirling.
  137.  * 
  138.  * Uses Lanczos-type approximation to ln(gamma) for z > 0. Reference: Lanczos,
  139.  * C. 'A precision approximation of the gamma function', J. SIAM Numer.
  140.  * Anal., B, 1, 86-96, 1964. Accuracy: About 14 significant digits except for
  141.  * small regions in the vicinity of 1 and 2.
  142.  * 
  143.  * Programmer: Alan Miller CSIRO Division of Mathematics & Statistics
  144.  * 
  145.  * Latest revision - 17 April 1988
  146.  * 
  147.  * Additions: Translated from fortran to C, code added to handle values z < 0.
  148.  * The global variable signgam contains the sign of the gamma function.
  149.  * 
  150.  * IMPORTANT: The signgam variable contains garbage until AFTER the call to
  151.  * lngamma().
  152.  * 
  153.  * Permission granted to distribute freely for non-commercial purposes only
  154.  * Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl
  155.  */
  156.  
  157. /* Local data, not visible outside this file 
  158.    static double   a[] =
  159.    {
  160.    0.9999999999995183E+00,
  161.    0.6765203681218835E+03,
  162.    -.1259139216722289E+04,
  163.    0.7713234287757674E+03,
  164.    -.1766150291498386E+03,
  165.    0.1250734324009056E+02,
  166.    -.1385710331296526E+00,
  167.    0.9934937113930748E-05,
  168.    0.1659470187408462E-06,
  169.    };   */
  170.  
  171. /* from Ray Toy */
  172. static double GPFAR a[] =
  173. {
  174.     .99999999999980993227684700473478296744476168282198,
  175.     676.52036812188509856700919044401903816411251975244084,
  176.     -1259.13921672240287047156078755282840836424300664868028,
  177.     771.32342877765307884865282588943070775227268469602500,
  178.     -176.61502916214059906584551353999392943274507608117860,
  179.     12.50734327868690481445893685327104972970563021816420,
  180.     -.13857109526572011689554706984971501358032683492780,
  181.     .00000998436957801957085956266828104544089848531228,
  182.     .00000015056327351493115583383579667028994545044040,
  183. };
  184.  
  185. static double lgamneg(z)
  186. double z;
  187. {
  188.     double tmp;
  189.  
  190.     /* Use reflection formula, then call lgampos() */
  191.     tmp = sin(z * PI);
  192.  
  193.     if (fabs(tmp) < MACHEPS) {
  194.     tmp = 0.0;
  195.     } else if (tmp < 0.0) {
  196.     tmp = -tmp;
  197.     signgam = -1;
  198.     }
  199.     return LNPI - lgampos(1.0 - z) - log(tmp);
  200.  
  201. }
  202.  
  203. static double lgampos(z)
  204. double z;
  205. {
  206.     double sum;
  207.     double tmp;
  208.     int i;
  209.  
  210.     sum = a[0];
  211.     for (i = 1, tmp = z; i < 9; i++) {
  212.     sum += a[i] / tmp;
  213.     tmp++;
  214.     }
  215.  
  216.     return log(sum) + LNSQRT2PI - z - 6.5 + (z - 0.5) * log(z + 6.5);
  217. }
  218.  
  219. static double lngamma(z)
  220. double z;
  221. {
  222.     signgam = 1;
  223.  
  224.     if (z <= 0.0)
  225.     return lgamneg(z);
  226.     else
  227.     return lgampos(z);
  228. }
  229.  
  230. # define GAMMA(x) lngamma ((x))
  231. #endif /* !GAMMA */
  232.  
  233. void f_erf()
  234. {
  235.     struct value a;
  236.     double x;
  237.  
  238.     x = real(pop(&a));
  239.  
  240. #ifdef HAVE_ERF
  241.     x = erf(x);
  242. #else
  243.     {
  244.     int fsign;
  245.     fsign = x >= 0 ? 1 : 0;
  246.     x = igamma(0.5, (x)*(x));
  247.     if (x == -1.0) {
  248.         undefined = TRUE;
  249.         x = 0.0;
  250.     } else {
  251.         if (fsign == 0)
  252.         x = -x;
  253.     }
  254.     }
  255. #endif
  256.     push(Gcomplex(&a, x, 0.0));
  257. }
  258.  
  259. void f_erfc()
  260. {
  261.     struct value a;
  262.     double x;
  263.  
  264.     x = real(pop(&a));
  265. #ifdef HAVE_ERFC
  266.     x = erfc(x);
  267. #else
  268.     {
  269.     int fsign;
  270.     fsign = x >= 0 ? 1 : 0;
  271.     x = igamma(0.5, (x)*(x));
  272.     if (x == 1.0) {
  273.         undefined = TRUE;
  274.         x = 0.0;
  275.     } else { 
  276.         x = fsign > 0 ? 1.0 - x : 1.0 + x ;
  277.     }
  278.     }
  279. #endif
  280.     push(Gcomplex(&a, x, 0.0));
  281. }
  282.  
  283. void f_ibeta()
  284. {
  285.     struct value a;
  286.     double x;
  287.     double arg1;
  288.     double arg2;
  289.  
  290.     x = real(pop(&a));
  291.     arg2 = real(pop(&a));
  292.     arg1 = real(pop(&a));
  293.  
  294.     x = ibeta(arg1, arg2, x);
  295.     if (x == -1.0) {
  296.     undefined = TRUE;
  297.     push(Ginteger(&a, 0));
  298.     } else
  299.     push(Gcomplex(&a, x, 0.0));
  300. }
  301.  
  302. void f_igamma()
  303. {
  304.     struct value a;
  305.     double x;
  306.     double arg1;
  307.  
  308.     x = real(pop(&a));
  309.     arg1 = real(pop(&a));
  310.  
  311.     x = igamma(arg1, x);
  312.     if (x == -1.0) {
  313.     undefined = TRUE;
  314.     push(Ginteger(&a, 0));
  315.     } else
  316.     push(Gcomplex(&a, x, 0.0));
  317. }
  318.  
  319. void f_gamma()
  320. {
  321.     register double y;
  322.     struct value a;
  323.  
  324.     y = GAMMA(real(pop(&a)));
  325.     if (y > 88.0) {
  326.     undefined = TRUE;
  327.     push(Ginteger(&a, 0));
  328.     } else
  329.     push(Gcomplex(&a, signgam * gp_exp(y), 0.0));
  330. }
  331.  
  332. void f_lgamma()
  333. {
  334.     struct value a;
  335.  
  336.     push(Gcomplex(&a, GAMMA(real(pop(&a))), 0.0));
  337. }
  338.  
  339. #ifndef BADRAND
  340.  
  341. void f_rand()
  342. {
  343.     struct value a;
  344.  
  345.     push(Gcomplex(&a, ranf(real(pop(&a))), 0.0));
  346. }
  347.  
  348. #else /* BADRAND */
  349.  
  350. /* Use only to observe the effect of a "bad" random number generator. */
  351. void f_rand()
  352. {
  353.     struct value a;
  354.  
  355.     static unsigned int y = 0;
  356.     unsigned int maxran = 1000;
  357.  
  358.     (void) real(pop(&a));
  359.     y = (781 * y + 387) % maxran;
  360.  
  361.     push(Gcomplex(&a, (double) y / maxran, 0.0));
  362. }
  363.  
  364. #endif /* BADRAND */
  365.  
  366. /** ibeta.c
  367.  *
  368.  *   DESCRIB   Approximate the incomplete beta function Ix(a, b).
  369.  *
  370.  *                           _
  371.  *                          |(a + b)     /x  (a-1)         (b-1)
  372.  *             Ix(a, b) = -_-------_--- * |  t     * (1 - t)     dt (a,b > 0)
  373.  *                        |(a) * |(b)   /0
  374.  *
  375.  *
  376.  *
  377.  *   CALL      p = ibeta(a, b, x)
  378.  *
  379.  *             double    a    > 0
  380.  *             double    b    > 0
  381.  *             double    x    [0, 1]
  382.  *
  383.  *   WARNING   none
  384.  *
  385.  *   RETURN    double    p    [0, 1]
  386.  *                            -1.0 on error condition
  387.  *
  388.  *   XREF      lngamma()
  389.  *
  390.  *   BUGS      none
  391.  *
  392.  *   REFERENCE The continued fraction expansion as given by
  393.  *             Abramowitz and Stegun (1964) is used.
  394.  *
  395.  * Permission granted to distribute freely for non-commercial purposes only
  396.  * Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl
  397.  */
  398.  
  399. static double ibeta(a, b, x)
  400. double a, b, x;
  401. {
  402.     /* Test for admissibility of arguments */
  403.     if (a <= 0.0 || b <= 0.0)
  404.     return -1.0;
  405.     if (x < 0.0 || x > 1.0)
  406.     return -1.0;;
  407.  
  408.     /* If x equals 0 or 1, return x as prob */
  409.     if (x == 0.0 || x == 1.0)
  410.     return x;
  411.  
  412.     /* Swap a, b if necessarry for more efficient evaluation */
  413.     return a < x * (a + b) ? 1.0 - confrac(b, a, 1.0 - x) : confrac(a, b, x);
  414. }
  415.  
  416. static double confrac(a, b, x)
  417. double a, b, x;
  418. {
  419.     double Alo = 0.0;
  420.     double Ahi;
  421.     double Aev;
  422.     double Aod;
  423.     double Blo = 1.0;
  424.     double Bhi = 1.0;
  425.     double Bod = 1.0;
  426.     double Bev = 1.0;
  427.     double f;
  428.     double fold;
  429.     double Apb = a + b;
  430.     double d;
  431.     int i;
  432.     int j;
  433.  
  434.     /* Set up continued fraction expansion evaluation. */
  435.     Ahi = gp_exp(GAMMA(Apb) + a * log(x) + b * log(1.0 - x) -
  436.          GAMMA(a + 1.0) - GAMMA(b));
  437.  
  438.     /*
  439.      * Continued fraction loop begins here. Evaluation continues until
  440.      * maximum iterations are exceeded, or convergence achieved.
  441.      */
  442.     for (i = 0, j = 1, f = Ahi; i <= ITMAX; i++, j++) {
  443.     d = a + j + i;
  444.     Aev = -(a + i) * (Apb + i) * x / d / (d - 1.0);
  445.     Aod = j * (b - j) * x / d / (d + 1.0);
  446.     Alo = Bev * Ahi + Aev * Alo;
  447.     Blo = Bev * Bhi + Aev * Blo;
  448.     Ahi = Bod * Alo + Aod * Ahi;
  449.     Bhi = Bod * Blo + Aod * Bhi;
  450.  
  451.     if (fabs(Bhi) < MACHEPS)
  452.         Bhi = 0.0;
  453.  
  454.     if (Bhi != 0.0) {
  455.         fold = f;
  456.         f = Ahi / Bhi;
  457.         if (fabs(f - fold) < fabs(f) * MACHEPS)
  458.         return f;
  459.     }
  460.     }
  461.  
  462.     return -1.0;
  463. }
  464.  
  465. /** igamma.c
  466.  *
  467.  *   DESCRIB   Approximate the incomplete gamma function P(a, x).
  468.  *
  469.  *                         1     /x  -t   (a-1)
  470.  *             P(a, x) = -_--- * |  e  * t     dt      (a > 0)
  471.  *                       |(a)   /0
  472.  *
  473.  *   CALL      p = igamma(a, x)
  474.  *
  475.  *             double    a    >  0
  476.  *             double    x    >= 0
  477.  *
  478.  *   WARNING   none
  479.  *
  480.  *   RETURN    double    p    [0, 1]
  481.  *                            -1.0 on error condition
  482.  *
  483.  *   XREF      lngamma()
  484.  *
  485.  *   BUGS      Values 0 <= x <= 1 may lead to inaccurate results.
  486.  *
  487.  *   REFERENCE ALGORITHM AS239  APPL. STATIST. (1988) VOL. 37, NO. 3
  488.  *
  489.  * Permission granted to distribute freely for non-commercial purposes only
  490.  * Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl
  491.  */
  492.  
  493. /* Global variables, not visible outside this file */
  494. static double pn1, pn2, pn3, pn4, pn5, pn6;
  495.  
  496. static double igamma(a, x)
  497. double a, x;
  498. {
  499.     double arg;
  500.     double aa;
  501.     double an;
  502.     double b;
  503.     int i;
  504.  
  505.     /* Check that we have valid values for a and x */
  506.     if (x < 0.0 || a <= 0.0)
  507.     return -1.0;
  508.  
  509.     /* Deal with special cases */
  510.     if (x == 0.0)
  511.     return 0.0;
  512.     if (x > XBIG)
  513.     return 1.0;
  514.  
  515.     /* Check value of factor arg */
  516.     arg = a * log(x) - x - GAMMA(a + 1.0);
  517.     if (arg < MINEXP)
  518.     return -1.0;
  519.     arg = gp_exp(arg);
  520.  
  521.     /* Choose infinite series or continued fraction. */
  522.  
  523.     if ((x > 1.0) && (x >= a + 2.0)) {
  524.     /* Use a continued fraction expansion */
  525.  
  526.     double rn;
  527.     double rnold;
  528.  
  529.     aa = 1.0 - a;
  530.     b = aa + x + 1.0;
  531.     pn1 = 1.0;
  532.     pn2 = x;
  533.     pn3 = x + 1.0;
  534.     pn4 = x * b;
  535.     rnold = pn3 / pn4;
  536.  
  537.     for (i = 1; i <= ITMAX; i++) {
  538.  
  539.         aa++;
  540.         b += 2.0;
  541.         an = aa * (double) i;
  542.  
  543.         pn5 = b * pn3 - an * pn1;
  544.         pn6 = b * pn4 - an * pn2;
  545.  
  546.         if (pn6 != 0.0) {
  547.  
  548.         rn = pn5 / pn6;
  549.         if (fabs(rnold - rn) <= GPMIN(MACHEPS, MACHEPS * rn))
  550.             return 1.0 - arg * rn * a;
  551.  
  552.         rnold = rn;
  553.         }
  554.         pn1 = pn3;
  555.         pn2 = pn4;
  556.         pn3 = pn5;
  557.         pn4 = pn6;
  558.  
  559.         /* Re-scale terms in continued fraction if terms are large */
  560.         if (fabs(pn5) >= OFLOW) {
  561.  
  562.         pn1 /= OFLOW;
  563.         pn2 /= OFLOW;
  564.         pn3 /= OFLOW;
  565.         pn4 /= OFLOW;
  566.         }
  567.     }
  568.     } else {
  569.     /* Use Pearson's series expansion. */
  570.  
  571.     for (i = 0, aa = a, an = b = 1.0; i <= ITMAX; i++) {
  572.  
  573.         aa++;
  574.         an *= x / aa;
  575.         b += an;
  576.         if (an < b * MACHEPS)
  577.         return arg * b;
  578.     }
  579.     }
  580.     return -1.0;
  581. }
  582.  
  583. /***********************************************************************
  584.      double ranf(double init)
  585.                 RANDom number generator as a Function
  586.      Returns a random floating point number from a uniform distribution
  587.      over 0 - 1 (endpoints of this interval are not returned) using a
  588.      large integer generator.
  589.      This is a transcription from Pascal to Fortran of routine
  590.      Uniform_01 from the paper
  591.      L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
  592.      with Splitting Facilities." ACM Transactions on Mathematical
  593.      Software, 17:98-111 (1991)
  594.  
  595.                GeNerate LarGe Integer
  596.      Returns a random integer following a uniform distribution over
  597.      (1, 2147483562) using the generator.
  598.      This is a transcription from Pascal to Fortran of routine
  599.      Random from the paper
  600.      L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
  601.      with Splitting Facilities." ACM Transactions on Mathematical
  602.      Software, 17:98-111 (1991)
  603. ***********************************************************************/
  604. static double ranf(init)
  605. double init;
  606. {
  607.  
  608.     long k, z;
  609.     static int firsttime = 1;
  610.     static long s1, s2;
  611.  
  612.     /* (Re)-Initialize seeds if necessary */
  613.     if (init < 0.0 || firsttime == 1) {
  614.     firsttime = 0;
  615.     s1 = 1234567890L;
  616.     s2 = 1234567890L;
  617.     }
  618.     /* Generate pseudo random integers */
  619.     k = s1 / 53668L;
  620.     s1 = Xa1 * (s1 - k * 53668L) - k * 12211;
  621.     if (s1 < 0)
  622.     s1 += Xm1;
  623.     k = s2 / 52774L;
  624.     s2 = Xa2 * (s2 - k * 52774L) - k * 3791;
  625.     if (s2 < 0)
  626.     s2 += Xm2;
  627.     z = s1 - s2;
  628.     if (z < 1)
  629.     z += (Xm1 - 1);
  630.  
  631.     /*
  632.      * 4.656613057E-10 is 1/Xm1.  Xm1 is set at the top of this file and is
  633.      * currently 2147483563. If Xm1 changes, change this also.
  634.      */
  635.     return (double) 4.656613057E-10 *z;
  636. }
  637.  
  638. /* ----------------------------------------------------------------
  639.    Following to specfun.c made by John Grosh (jgrosh@arl.mil)
  640.    on 28 OCT 1992.
  641.    ---------------------------------------------------------------- */
  642.  
  643. void f_normal()
  644. {                /* Normal or Gaussian Probability Function */
  645.     struct value a;
  646.     double x;
  647.  
  648.     /* ref. Abramowitz and Stegun 1964, "Handbook of Mathematical 
  649.        Functions", Applied Mathematics Series, vol 55,
  650.        Chapter 26, page 934, Eqn. 26.2.29 and Jos van der Woude 
  651.        code found above */
  652.  
  653.     x = real(pop(&a));
  654.  
  655.     x = 0.5 * SQRT_TWO * x;
  656. #ifdef HAVE_ERF
  657.     x = 0.5 * (1.0 + erf(x));
  658. #else
  659.     {
  660.     int fsign;
  661.     fsign = x >= 0 ? 1 : 0;
  662.     x = igamma(0.5, (x)*(x));
  663.     if (x == 1.0) {
  664.         undefined = TRUE;
  665.         x = 0.0;
  666.     } else { 
  667.         if (fsign == 0)
  668.         x = -(x);
  669.         x = 0.5 * (1.0 + x);
  670.     }
  671.     }
  672. #endif
  673.     push(Gcomplex(&a, x, 0.0));
  674. }
  675.  
  676. void f_inverse_normal()
  677. {                /* Inverse normal distribution function */
  678.     struct value a;
  679.     double x;
  680.  
  681.     x = real(pop(&a));
  682.  
  683.     if (x <= 0.0 || x >= 1.0) {
  684.     undefined = TRUE;
  685.     push(Gcomplex(&a, 0.0, 0.0));
  686.     } else {
  687.     push(Gcomplex(&a, inverse_normal_func(x), 0.0));
  688.     }
  689. }
  690.  
  691.  
  692. void f_inverse_erf()
  693. {                /* Inverse error function */
  694.     struct value a;
  695.     double x;
  696.  
  697.     x = real(pop(&a));
  698.  
  699.     if (fabs(x) >= 1.0) {
  700.     undefined = TRUE;
  701.     push(Gcomplex(&a, 0.0, 0.0));
  702.     } else {
  703.     push(Gcomplex(&a, inverse_error_func(x), 0.0));
  704.     }
  705. }
  706.  
  707. static double inverse_normal_func(p)
  708. double p;
  709. {
  710.     /* 
  711.        Source: This routine was derived (using f2c) from the 
  712.        FORTRAN subroutine MDNRIS found in 
  713.        ACM Algorithm 602 obtained from netlib.
  714.  
  715.        MDNRIS code contains the 1978 Copyright 
  716.        by IMSL, INC. .  Since MDNRIS has been 
  717.        submitted to netlib it may be used with 
  718.        the restriction that it may only be 
  719.        used for noncommercial purposes and that
  720.        IMSL be acknowledged as the copyright-holder
  721.        of the code.
  722.      */
  723.  
  724.     /* Initialized data */
  725.     static double eps = 1e-10;
  726.     static double g0 = 1.851159e-4;
  727.     static double g1 = -.002028152;
  728.     static double g2 = -.1498384;
  729.     static double g3 = .01078639;
  730.     static double h0 = .09952975;
  731.     static double h1 = .5211733;
  732.     static double h2 = -.06888301;
  733.     static double sqrt2 = 1.414213562373095;
  734.  
  735.     /* Local variables */
  736.     static double a, w, x;
  737.     static double sd, wi, sn, y;
  738.  
  739.     /* Note: 0.0 < p < 1.0 */
  740.  
  741.     /* p too small, compute y directly */
  742.     if (p <= eps) {
  743.     a = p + p;
  744.     w = sqrt(-(double) log(a + (a - a * a)));
  745.  
  746.     /* use a rational function in 1.0 / w */
  747.     wi = 1.0 / w;
  748.     sn = ((g3 * wi + g2) * wi + g1) * wi;
  749.     sd = ((wi + h2) * wi + h1) * wi + h0;
  750.     y = w + w * (g0 + sn / sd);
  751.     y = -y * sqrt2;
  752.     } else {
  753.     x = 1.0 - (p + p);
  754.     y = inverse_error_func(x);
  755.     y = -sqrt2 * y;
  756.     }
  757.     return (y);
  758. }
  759.  
  760.  
  761. static double inverse_error_func(p)
  762. double p;
  763. {
  764.     /* 
  765.        Source: This routine was derived (using f2c) from the 
  766.        FORTRAN subroutine MERFI found in 
  767.        ACM Algorithm 602 obtained from netlib.
  768.  
  769.        MDNRIS code contains the 1978 Copyright 
  770.        by IMSL, INC. .  Since MERFI has been 
  771.        submitted to netlib, it may be used with 
  772.        the restriction that it may only be 
  773.        used for noncommercial purposes and that
  774.        IMSL be acknowledged as the copyright-holder
  775.        of the code.
  776.      */
  777.  
  778.  
  779.  
  780.     /* Initialized data */
  781.     static double a1 = -.5751703;
  782.     static double a2 = -1.896513;
  783.     static double a3 = -.05496261;
  784.     static double b0 = -.113773;
  785.     static double b1 = -3.293474;
  786.     static double b2 = -2.374996;
  787.     static double b3 = -1.187515;
  788.     static double c0 = -.1146666;
  789.     static double c1 = -.1314774;
  790.     static double c2 = -.2368201;
  791.     static double c3 = .05073975;
  792.     static double d0 = -44.27977;
  793.     static double d1 = 21.98546;
  794.     static double d2 = -7.586103;
  795.     static double e0 = -.05668422;
  796.     static double e1 = .3937021;
  797.     static double e2 = -.3166501;
  798.     static double e3 = .06208963;
  799.     static double f0 = -6.266786;
  800.     static double f1 = 4.666263;
  801.     static double f2 = -2.962883;
  802.     static double g0 = 1.851159e-4;
  803.     static double g1 = -.002028152;
  804.     static double g2 = -.1498384;
  805.     static double g3 = .01078639;
  806.     static double h0 = .09952975;
  807.     static double h1 = .5211733;
  808.     static double h2 = -.06888301;
  809.  
  810.     /* Local variables */
  811.     static double a, b, f, w, x, y, z, sigma, z2, sd, wi, sn;
  812.  
  813.     x = p;
  814.  
  815.     /* determine sign of x */
  816.     if (x > 0)
  817.     sigma = 1.0;
  818.     else
  819.     sigma = -1.0;
  820.  
  821.     /* Note: -1.0 < x < 1.0 */
  822.  
  823.     z = fabs(x);
  824.  
  825.     /* z between 0.0 and 0.85, approx. f by a 
  826.        rational function in z  */
  827.  
  828.     if (z <= 0.85) {
  829.     z2 = z * z;
  830.     f = z + z * (b0 + a1 * z2 / (b1 + z2 + a2
  831.                      / (b2 + z2 + a3 / (b3 + z2))));
  832.  
  833.     /* z greater than 0.85 */
  834.     } else {
  835.     a = 1.0 - z;
  836.     b = z;
  837.  
  838.     /* reduced argument is in (0.85,1.0), 
  839.        obtain the transformed variable */
  840.  
  841.     w = sqrt(-(double) log(a + a * b));
  842.  
  843.     /* w greater than 4.0, approx. f by a 
  844.        rational function in 1.0 / w */
  845.  
  846.     if (w >= 4.0) {
  847.         wi = 1.0 / w;
  848.         sn = ((g3 * wi + g2) * wi + g1) * wi;
  849.         sd = ((wi + h2) * wi + h1) * wi + h0;
  850.         f = w + w * (g0 + sn / sd);
  851.  
  852.         /* w between 2.5 and 4.0, approx. 
  853.            f by a rational function in w */
  854.  
  855.     } else if (w < 4.0 && w > 2.5) {
  856.         sn = ((e3 * w + e2) * w + e1) * w;
  857.         sd = ((w + f2) * w + f1) * w + f0;
  858.         f = w + w * (e0 + sn / sd);
  859.  
  860.         /* w between 1.13222 and 2.5, approx. f by 
  861.            a rational function in w */
  862.     } else if (w <= 2.5 && w > 1.13222) {
  863.         sn = ((c3 * w + c2) * w + c1) * w;
  864.         sd = ((w + d2) * w + d1) * w + d0;
  865.         f = w + w * (c0 + sn / sd);
  866.     }
  867.     }
  868.     y = sigma * f;
  869.     return (y);
  870. }
  871.